home *** CD-ROM | disk | FTP | other *** search
/ The World of Computer Software / The World of Computer Software.iso / yamp16.zip / TESTREG.CPP < prev    next >
C/C++ Source or Header  |  1993-01-11  |  6KB  |  264 lines
  1.  
  2. /*
  3. YAMP - Yet Another Matrix Program
  4. Version: 1.6
  5. Author: Mark Von Tress, Ph.D.
  6. Date: 01/11/93  
  7.  
  8. Copyright(c) Mark Von Tress 1993
  9. Portions of this code are (c) 1991 by Allen I. Holub and are used by
  10. permission
  11.  
  12. DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
  13. WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
  14. TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
  15. ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
  16. FROM USE OF THIS PROGRAM.
  17.  
  18. */
  19.  
  20. #include "virt.h"
  21.  
  22. //required global declaration for the
  23. //  matrix stack object
  24.  
  25. unsigned int _stklen = STACKLENGTH;
  26. MStack *Dispatch = new MStack;
  27.  
  28. VMatrix& ExactHilbInv( int n )
  29. {
  30.    // construct exact hilbert matrix inverse
  31.    Dispatch->Inclevel();
  32.    VMatrix Hi("Hilbert Inverse",n,n);
  33.  
  34.    double diag, dn = n;
  35.    for( int i=1; i<=n; i++){
  36.       double im1 = (double) (i-1.0);
  37.       diag = ( i == 1 ) ? dn : ((dn-im1)*diag*(dn+im1))/(im1*im1);
  38.  
  39.       long double rest = diag*diag;
  40.       Hi.M(i,i) = rest/(2.0*im1+1.0);
  41.       for( int j=i+1; j<=n; j++ ){
  42.          double jm1 = (double) (j-1);
  43.          rest = -((dn-jm1)*rest*(dn+jm1))/(jm1*jm1);
  44.          Hi.M(i,j) = rest/(im1+jm1+1.0);
  45.          Hi.M(j,i) = Hi.m(i,j);
  46.       }
  47.    }
  48.    Hi.DisplayMat();
  49.    Dispatch->Push( Hi );
  50.    return Dispatch->DecReturn();
  51. }
  52.  
  53.  
  54. void TestHilb( int hilb_ind )
  55. // test svd on hilbert matrix
  56. {
  57.   Dispatch->Inclevel();
  58.   VMatrix hilb = Kron( Fill(1,hilb_ind,1), Index(1, hilb_ind));
  59.   // h(i,j) = 1.0/(i+j-1)
  60.   hilb = 1.0 / ( (hilb+Tran(hilb))-1.0);
  61.   hilb.Nameit("Hilbert Matrix");
  62.   hilb.DisplayMat();
  63.   // use formula. Looses accuracy of 0.001 at n>=11
  64.   (ExactHilbInv( hilb_ind )*hilb ).DisplayMat();
  65.  
  66.   // use sweep
  67.   VMatrix Gauss = Inv(hilb)*hilb;
  68.   Gauss.Nameit("Gaussian elimination");
  69.   Gauss.DisplayMat();
  70.  
  71.   // use singular valued decomposition
  72.   VMatrix S, V, D;
  73.   Svd( hilb, S, V, D);
  74.   VMatrix t = Fill(V.r, V.r, 0.0);
  75.   for (int i = 1; i <= V.r; i++) {
  76.     double tt = V.m(i, 1);
  77.     t.M(i, i) = 1.0 / tt;
  78.   }
  79.  
  80.   VMatrix Ggauss = (D*t*Tran(S))*hilb;
  81.   Ggauss.Nameit("Generalized inversion");
  82.   Ggauss.DisplayMat();
  83.  
  84.   Dispatch->Declevel();
  85. }
  86.  
  87. ///////////////////// now test regressions on sample data
  88.  
  89. VMatrix &getx(int N)
  90. // create an x matrix
  91. {
  92.     Dispatch->Inclevel();
  93.  
  94.     VMatrix x = Index( 1, N ) - N/2;
  95.     VMatrix c1 = Fill(N, 1, 1.0);
  96.     VMatrix x2 = x % x;
  97.     x = Ch(c1, Ch(x, x2));
  98.  
  99.     // push x onto the stack
  100.     Dispatch->Push(x);
  101.     // decrement the subroutine nesting level
  102.     // and return the stack top
  103.     return Dispatch->DecReturn();
  104. }
  105. #ifndef random
  106. #define random( max ) ((rand() % (int)(((max)+1) - (0))) + (0))
  107. #endif
  108.  
  109. VMatrix &gety(VMatrix &x, VMatrix &beta)
  110. // create a y matrix
  111. {
  112.     Dispatch->Inclevel();
  113.  
  114.     VMatrix y = x*beta;
  115.     srand(123);
  116.     for (int i = 1; i <= y.r; i++) {
  117.         // use sum of 3 uniforms for an approximate
  118.         // normal random variable
  119.         int u = random(100) + random(100) + random(100) + 3;
  120.         y.M(i, 1) = y.m(i, 1) + ((double) (u - 150)) / 300.0;
  121.     }
  122.     Dispatch->Push(y);
  123.     return Dispatch->DecReturn();
  124. }
  125.  
  126. VMatrix ®ression(VMatrix& x, VMatrix& y)
  127. // do a multiple linear regression
  128. {
  129.     Dispatch->Inclevel();
  130.  
  131.     int N = x.r, p = x.c;
  132.     // solve for regression parameters using sweep
  133.     VMatrix yx = Ch(y, x);
  134.     VMatrix reg = Sweep(2, p + 1, Tran(yx) *yx);
  135.     // calculate mean square error
  136.     reg.M(1, 1) = reg.m(1, 1) / ((double) (N - p));
  137.     reg.DisplayMat();
  138.  
  139.  
  140.     // solve regression using normal equations
  141.     VMatrix betahat = Inv(Tran(x) *x) *Tran(x) *y;
  142.  
  143.     Dispatch->Push(betahat);
  144.     return Dispatch->DecReturn();
  145. }
  146. VMatrix &QRregression(VMatrix &x, VMatrix &y)
  147. {
  148.     // use QR decomposition to solve regression
  149.  
  150.     Dispatch->Inclevel();
  151.     VMatrix Q, R;
  152.  
  153.     QR(x, Q, R);
  154.     VMatrix betahat = Inv(Tran(R) *R) *Tran(R) *Tran(Q) *y;
  155.  
  156.     Dispatch->Push(betahat);
  157.     return Dispatch->DecReturn();
  158. }
  159.  
  160. VMatrix &GinvRegression(VMatrix &x, VMatrix &y)
  161. {
  162.     // use Ginv to solve normal equations
  163.     Dispatch->Inclevel();
  164.  
  165.     VMatrix betahat = Ginv(Tran(x) *x) *Tran(x) *y;
  166.  
  167.     Dispatch->Push(betahat);
  168.     return Dispatch->DecReturn();
  169. }
  170.  
  171.  
  172. VMatrix &SVDregression(VMatrix &x, VMatrix &y)
  173. {
  174.     // use SVD to solve regression x = S Diag(V) Tran(D)
  175.     Dispatch->Inclevel();
  176.     VMatrix S, V, D;
  177.  
  178.     Svd(x, S, V, D);
  179.     VMatrix t = Fill(V.r, V.r, 0.0);
  180.     for (int i = 1; i <= V.r; i++) {
  181.         double tt = V.m(i, 1);
  182.         t.M(i, i) = (fabs(tt) <= 0.0) ? 0.0 : 1.0 / tt;
  183.     }
  184.     VMatrix betahat = D*t*Tran(S)*y;
  185.  
  186.     Dispatch->Push(betahat);
  187.     return Dispatch->DecReturn();
  188. }
  189.  
  190. VMatrix &GetSerialCovar( VMatrix &R )
  191. {
  192.    // Parameters to CORR in Timeslab
  193.    // correlogram = CORR(x=R,n=R.r,M=R.r-1,Q=2*R.r,
  194.    //                    iopt=1,r0=r0, per = spectdens)
  195.    // covar = correlogram*r0
  196.    Dispatch->Inclevel();
  197.    VMatrix centered, z, covar, spectdens;
  198.    double n = (double) R.r;
  199.    // center a column vector
  200.    centered = R - Sum(R).m(1, 1) / n;
  201.    // zero pad to length 2n: 2n periodic for full
  202.    // sample spectral density
  203.    centered = Cv(centered, Fill(R.r, R.c, 0));
  204.    // take fft
  205.    z = Fft(centered);
  206.    // take convolution : gives sample spectral density
  207.    spectdens = Sum(z % z / n, COLUMNS);
  208.    // inverse fft for serial correlation (autocorrelation)
  209.    covar = Fft(spectdens, FFALSE);
  210.    // throw away last half.
  211.    covar = Submat(covar, R.r, 2);
  212.    Dispatch->Push(covar);
  213.    return Dispatch->DecReturn();
  214. }
  215.  
  216.  
  217. main()
  218. {
  219.     Dispatch->Inclevel();
  220.     VMatrix x, beta("beta", 3, 1), y, betahat, resids, serial;
  221.  
  222.     TestHilb( 11 );
  223.  
  224.     Setwid(15);
  225.     Setdec(10);
  226.  
  227.  
  228.     beta.M(1, 1) = 1;
  229.     beta.M(2, 1) = 0.5;
  230.     beta.M(3, 1) = 0.25;
  231.  
  232.     x = getx(50);
  233.     y = gety(x, beta);
  234.  
  235.     betahat = regression(x, y);
  236.     betahat.Nameit("Text book betahat");
  237.     (Tran(beta)).DisplayMat();
  238.     (Tran(betahat)).DisplayMat();
  239.  
  240.     betahat = QRregression(x, y);
  241.     betahat.Nameit("QR betahat");
  242.     (Tran(beta)).DisplayMat();
  243.     (Tran(betahat)).DisplayMat();
  244.  
  245.     betahat = GinvRegression(x, y);
  246.     betahat.Nameit("Ginv betahat");
  247.     (Tran(beta)).DisplayMat();
  248.     (Tran(betahat)).DisplayMat();
  249.  
  250.     betahat = SVDregression(x, y);
  251.     betahat.Nameit("SVD regression");
  252.     (Tran(beta)).DisplayMat();
  253.     (Tran(betahat)).DisplayMat();
  254.  
  255.     resids = y - x*betahat;
  256.     serial = GetSerialCovar(resids);
  257.     (Tran(Submat(serial, 5, 1))).DisplayMat();
  258.  
  259.     Setwid(6);
  260.     Setdec(3);
  261.  
  262.     vclose();
  263. }
  264.